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Using T — Monte Carlo simulation, we study the relax- 
ation of graph coloring (K-COL) and satisfiability (K-SAT), 
two hard problems that have recently been shown to possess 
a phase transition in solvability as a parameter is varied. A 
change from exponentially fast to power law relaxation, and a 
transition to freezing behavior are found. These changes take 
place for smaller values of the parameter than the solvability 
transition. Results for the coloring problem for colorable and 
clustered graphs and for the fraction of persistent spins for 
satisfiability are also presented. 

PACS numbers: 75.10.Hk (Classical spin models), 75.10.Nr 
(Spin-glass and other random models), 02.10.Eb (Combina- 
torics), 89.80.+h (Computer science and technology) 



I. INTRODUCTION 

Computers have made it possible for physicists to do 
experiments without leaving their offices. By simulating 
systems from Nature, or simple models from theoreti- 
cal physics, a tremendous amount of information can be 
gained. The advance of computer hardware and software 
has made it possible to simulate systems consisting of 
vast numbers of interacting particles. Physics, in turn, 
has inspired new algorithms for solving problems in com- 
puter science. For instance, the method of simulated an- 
nealing Q, based on Monte Carlo simulations with the 
Metropolis algorithm j|, can sometimes find good solu- 
tions much quicker than traditional algorithms S. 

Physicists have also started to study systems that are 
not found in Nature, but instead come from computer 
science 0^|. As an example, phase transitions in op- 
timization problems have been discovered and studied 
using statistical mechanics ||J^] . Other problems studied 
using physical methods include the knapsack-problem || , 
graph-partitioning Q , minimax-gamcs f|3~0Tl , the 8-Queens 
problem [0, number partitioning fl2] , |13H and the stable 
marriage problem Field theory has also been used 
to study, e.g., the enumeration of Hamiltonian cycles on 
graphs pi and coloring of random, planar graphs p6| , p^7t . 

Here we present a study of the T = relaxation behav- 
ior of hard optimization problems. Using Monte Carlo 
simulation, we have measured the energy of a system 



starting in a random (excited) state and slowly relaxing 
into the ground state or a low-lying excited state. We 
find qualitatively different behavior for hard and easy in- 
stances of the problems. One of many models (e.g., [^8|) 
for which relaxation has been studied is the ferromagnetic 
Ising model on a regular lattice. The models studied here 
differ from the Ising model in several respects. They are 
random, i.e., there is no pattern in the interactions be- 
tween the spins, and the interactions have infinite range. 
K-COL can be viewed as a Potts model on a random 
graph with finite connectivity. Another important differ- 
ence is that these models can suffer from frustration. 

This article is organized in the following way: Sec- 
tion [n] introduces some concepts from computer science, 
while the problems we study are described in section IH . 
Previous work on the solvability transition is reviewed 
in section IV and some approximate explanations in sec- 
tion [v|. Section VI describes our results on the relaxation 



behavior, section VH compares the results to other local 
search methods. Additional measurements, e.g., of the 
fraction of persistent spins, are found in section VIII. 



Section IX discusses possible explanations for the relax- 
ation behavior and comments on the importance of en- 
tropy barriers. Conclusions and a discussion are con- 
tained in section |x|. 



II. COMPUTER SCIENCE FOR PHYSICISTS 

Computer scientists classify problems according to the 
maximal amount of resources needed for their solution. 
The most important resource is time, but it is also pos- 
sible to distinguish between problems that require qual- 
itatively different amounts of memory. For example, a 
list of N elements can always be sorted in time less 
than kN log N, where k is some constant ]l9J . The 
problems whose running time on a universal Turing ma- 
chine (e.g., pp[|) is bounded by a polynomial in their 
size are said to be in the class P. The important class 
NP (for non-deterministic polynomial) consists of those 
problems where it can be checked in polynomial time 
whether a proposed solution actually solves the problem. 
(A non- deterministic Turing machine would be able to 
solve NP problems in polynomial time.) It is obvious 
that P C NP, but there is no proof that P ^ NP. 
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However, most people believe that there are NP prob- 
lems whose worst-case instances take exponential time to 
solve on a universal Turing machine. 

The class NP-complctc (or NPC) arc the most impor- 
tant problems in NP. A problem of size N is in NPC 
if all other NP problems can be transformed into it in 
time at most polynomial in N. A method to solve an 
NPC problem efficiently can thus be used to solve any 
NP problem efficiently. It is known that if P ^ NP 
then there are problems in NP that are in neither P nor 
NPC. A problem is called NP-hard if it is at least as 
difficult as the most difficult NP problems; NPC is the 
intersection of NP and NP-hard. 

It is worth emphasizing that it is the worst case com- 
plexity that determines whether a problem is in NPC. 
The average time needed to solve an NPC problem 
(given a distribution of problems) may still be polyno- 
mial in problem size. The properties of K-COL and K- 
SAT studied here are related to average case behavior, 
and do not address the question of whether P ^ NP. 

A nice introduction to most elementary concepts from 
theoretical computer science can be found in ^T|. A 
modern reference on complexity theory and NP problems 
is p3), while p3| has an extensive list of NPC problems. 



two clauses and three variables. Applications of K-SAT 
include theorem proving, VLSI design, and learning. 

In K-SAT, each clause forbids one of the 2 K possible 
assignments for its variables. In the same way, an edge 
in a graph forbids K of the K 2 different colorings of its 
vertices. For both problems, there are M constraints 
on the solutions. The energy e of a problem instance 
is defined as the number of unsatisfied constraints per 
variable. Both K-SAT and K-COL are in P for K = 
2 and in NPC for K > 3 |§f]. The related problem 
(MAX-K-SAT) of trying to minimize the number of 
unsatisfied clauses in K-SAT is in NPC even for K = 2. 

The most important problems are those where the 
number of constraints is of the same order as the number 
of variables, M = aN. The scheduling problem described 
above fulfills this condition, for example. For graph col- 
oring, a — 7/2, where 7 is the connectivity of the graph. 
The connectivity (or average degree) is defined as the 
mean number of edges exiting each node. For a graph 
with n vertices and e edges, it is 2e/n. 

Below, we will use a for K-SAT and 7 when we talk of 
K-COL. We will concentrate on K-COL; most results 
for K-SAT are similar. 



III. K-SAT AND K-COL 

Two important problems in NPC are graph coloring 
(K-COL) and satisfiability testing (SAT). Graph color- 
ing is the problem of coloring a graph with N vertices 
and M edges using K colors so that no two adjacent ver- 
tices have the same color. In physical terms, K-COL is 
the problem of finding a ground state without frustrated 
bonds in an antiferromagnetic if-state Potts model on 
a random graph. Related models have been studied by 
Baillie, Johnston, and coworkers (e.g. p4| , p5| ]), who con- 
sidered Potts models on 0™-model Feynman diagrams. 
They found similarities between models on (j) 3 and (ft 
graphs and Bethe lattices, and showed that mean-field 
theories work well for describing both ferromagnetic and 
antiferromagnetic models on Feynman diagrams. 

The most natural application of graph coloring is in 
scheduling problems. For example, a school where each 
teacher and student can be involved in several different 
classes must schedule the classes so that no collisions oc- 
cur. If there are K different time slots available, this is 
K-COL. 

Satisfiability was the first problem shown to be in 
NPC ]2^1 - It is the problem of finding an assignment 
of true or false to N variables so that a boolean formula 
in them is satisfied. In K-SAT, this formula is written 
in conjunctive normal form (CNF), that is, it consists 
of the logical AND of M clauses, each clause being the 
OR of K (possibly negated) variables, where the same 
clause may appear more than once in a formula. For ex- 
ample, [x V y) A (y V -iz) is an instance of 2-SAT with 



IV. THE TRANSITION 

For physicists, an interesting property of these prob- 
lems is that they contain phase transitions |Q,^,^8). As 
the number of constraints increases, there is a transition 
from a region where almost all instances of the problem 
are solvable to a region where practically none can be 
solved. Physically, the transition is from a region where 
the ground state has zero energy to one where it is finite. 
For K = 2 the transition can be seen as a transition 
between a P problem (finding a perfect solution to the 
problem) and an NPC problem (finding an assignment 
of variables that minimizes the number of unsatisfied con- 
straints). 

An approximation similar to mean-field theory has 
been used by Williams and Hogg (e.g., p9|| ) and others 
to explain some of the properties of the phase transition. 
Recently, Friedgut has also made some progress towards 
showing rigorously the existence of a sharp transition in 
solvability for both K-SAT §| and K-COL ||. 

Related to this phase transition in problem solvability, 
there is a transition in how difficult it is to solve a prob- 
lem or show that no solutions exist |27 3^j3^]. If there 
are few constraints on the solution (the problem is under- 
constrained), it is easy to find one. Similarly, if there are 
so many constraints that the problem is overconstrained, 
not much effort is needed to show that it is unsolvable. In 
between these regions, where the problems are critically 
constrained, there is a peak in problem difficulty. This is 
called the "easy- hard-easy" transition p2| . 

K-SAT has recently been studied by a number of 
physicists. Kirkpatrick and Selman studied the phase 
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transition using finite-size scaling methods, and Monas- 
son and Zecchina H| used the replica method to 
show that the entropy of K-SAT stays finite at the tran- 
sition. This means that below the transition there are 
several solutions to each problem, all of which develop 
inconsistencies as the critical a is passed. Another prob- 
lem that has been studied using statistical mechanics is 
the number partitioning problem. Mertens has recently 
shown that it too has a phase transition ]l3] ]. The rele- 
vant parameter here is the ratio between the number of 
bits of input data and the number of variables. 

Another problem in NPC that also shows a transi- 
tion |3(| is the traveling salesperson problem (TSP), 
where the objective is to find a tour of minimum length 
visiting N given distinct cities. A difficulty in studying 
this problem is that there is no natural parameter (like 
a and 7) that distinguishes between under- and over- 
constrained problems. To get one, the TSP must be 
reformulated as a decision problem: is there a Hamilto- 
nian path of length less than 11 The parameter I plays 
the same role as a — for a given distribution of problems 
there is an l c such that HI ^> l c , almost all instances have 
a tour with length < I, but if I <^ l c practically no such 
tours exist. Traditionally, most NPC problems are for- 
mulated as decision rather than optimization problems. 

There are many NPC problems that contain no ob- 
vious parameter which makes it difficult to say if the 
solvability phase transition exists in all NPC-problcms 
or in just a few. There have been attempts to formulate 
a more general parameter (e.g., [p7|l ), with the drawback 
that it requires us to approximate the number of solu- 
tions. 

Phase transitions have also been found in problems 
beyond NPC, e.g., in QSAT |3q| , a harder version of 
satisfiability where the boolean variables are quantified 
by either V or 3 (in ordinary SAT, all variables are 
existentially quantified). This problem is known to be 
PSPACE-completc f2^| , meaning that it is at least as 
hard as all problems that can be solved by a universal 
Turing machine without time limits but using memory 
at most polynomial in problem size. 



V. AN APPROXIMATE THEORY FOR THE 
TRANSITION 

The approximation proposed by Williams and 
Hogg [||J2^] assumes that the constraints in the problem 
are independent. In physical terms it simply means per- 
forming an annealed rather than a quenched average over 
the disorder. It is exact for graphs without loops and for 
satisfiability problems where no variable is contained in 
more than one clause. The probability that an indepen- 
dent constraint is violated is p = 1/2 K for K-SAT and 
p = K/K 2 = 1/K for K-COL. The probability of hav- 
ing none of M constraints violated can be approximated 
as (1 — p) M , ignoring correlations between constraints, 



such as triangles in graphs. The number of solutions for 
K-COL is then 
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N soi = K<\i--y N '\ 

K 

and for K-SAT the expression is 

iv sol - 2 N (i - -^r N . 



(i) 



(2) 



Using the inclusion-exclusion principle it is possible to 
write an exact expression for jV so i |6|. The inclusion- 
exclusion principle is the generalization of the simple for- 
mula 

P(A U B) = P(A) + P{B) - P(A n B) 

from mathematical statistics. If we let Ai be the event 
that constraint i is violated, it expresses the probability 
that any (i.e., at least one) constraint is violated in terms 
of the probabilities of one, two, three or more constraints 
being violated simultaneously 



P(U i A l )=J2(-l) r+1 Sr 



(3) 



where S r is the probability of exactly r constraints being 
violated simultaneously. The number of solutions can 
now be found as 



Wad = iV ta t(l--P(Ui^)) > 



(4) 



where AT to t is the number of possible assignments of the 
variables, N tot = K N for K-COL and N tot = 2 N for 
K-SAT. For K-COL, Si = MR' 1 , since there are M 
edges and each of them eliminates K N_1 (of the K ) 
solutions. For S2, we need to express the number of states 
that are eliminated by each of two edges. This is given 

by ( ) K~ 2 , while the expression for 



S 3 = 



M 
3 



K 
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{K~ z - K )t, 



(5) 



requires knowledge of the number of triangles, t, in the 
graph. Expression (||) can be understood by noting that 
if two edges in a triangle are frustrated, the third is al- 
ways frustrated too. It can be shown that t is Poisson- 
distributed with mean 7 3 /6. To calculate Si for i > 4, 
we also need to know the distribution of more complex 
sub-graphs. 

The critical value of the parameter can be approxi- 
mated as that 7 which gives iV so i = 1 in (|l]), giving 



7c = -2 



log if 



iog(i-£; 



(6) 



For K = 3, equation gives 7 C = 5.4 for K-COL 
and a c = 5.2 for K-SAT. These values are larger than 
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the experimental values of j c = 4.6 and a c — 4.17. For 
K-SAT, this approximation has been independently in- 
troduced several times |39|, and references therein] . 

This calculation of the critical value of 7 ignores all cor- 
relations between different constraints in the problem. It 
gives an upper bound for 7 C and is analogous to studying 
a forest, a graph without cycles, in which all edges are 
violated with a probability p. Taking correlations into 
account reduces the number of solutions ||] . 

Kirousis et al |39| have introduced a new method to get 
an upper bound for the critical parameter of K-SAT. 
For K = 3, they prove a c < 4.598, and it is principle 
possible to get better bounds by including more terms in 
their expansion. A lower bound has also been found jfO), 
a c > 3.003. It is however difficult to generalize these 
methods to other problems, such as K-COL. 



VI. THE RELAXATION BEHAVIOR 

We have studied the relaxation of the energy e, de- 
fined as the number of unsatisfied constraints per spin, 
of K-SAT and K-COL using T = Monte Carlo (MC) 
simulations and the Metropolis single-flip algorithm 0. 
A simple case where the relaxation can be understood is 
the ferromagnetic Ising model on a regular lattice. For 
this model, the energy decreases as e ~ t -1 / 3 if the order 
parameter is conserved by the dynamics, while e ~ i -1 / 2 
if single spin flips that allow the total magnetization to 
change are used. These forms of relaxation behavior can 
be explained by noting that spins with the same orien- 
tation will cluster and form domains (e.g., [Q) with a 
well-defined energy. This explanation does not immedi- 
ately carry over to our problems, since there is no known 
simple expression for the energy as a function of a length 
scale. 

For the NPC problems studied here, the energy can 
not be expected always to reach zero, since there are no 
perfect solutions to problems with many constraints. The 
T = MC algorithm can also get stuck in local minima. 
A small graph where this can happen for 3-COL is shown 
in figure [l|. 
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this leaves the energy unchanged or lowers it, otherwise 
the spin is left unchanged. In temperature T > sim- 
ulations, a flip that raises the energy A units is allowed 
with probability exp[— A/T]. 

For K-COL, we generated the problems by randomly 

' N 



selecting M distinct edges from the 
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N vertices. (In graph terminology |4]f] , this corresponds 
to using the Q(N, M) model for random graphs.) For K- 
SAT, each clause was generated by selecting K variables, 
where each variable was negated with probability \ . The 
formula was then generated by performing this process 
M times. Clauses with repeated variables were allowed 
in the expressions, and also repetition of clauses. 

In all random systems, the question of whether or not 
the measured quantities are self-averaging is important. 
Self-averaging means that the properties of the entire (in- 
finite) system can be understood in terms of the proper- 
ties of its local subsystems. For random systems with lo- 
cal interactions, there is a simple argument for this (see 
e.g., 0)), but if the interactions are global, the situ- 
ation is more complex. However, for spin glasses such 
as the Sherrington-Kirkpatrick model, the energy and 
other simple quantities are self-averaging, and we as- 
sume that the energy is self-averaging also in K-SAT 
and K-COL. This assumption of self-averaging is sup- 
ported by our simulations; averaging over a small number 
of large graphs appears to give the same results as averag- 
ing over many small graphs. Schreiber and Martin |43] 
have recently studied various local search methods for 
the graph partitioning problem and found strong numer- 
ical evidence for self-averaging even for sparse graphs. 
They also provide some arguments for why self-averaging 
should hold, and conjecture that their result holds for all 
constraint satisfaction problems. 

Our simulations show a transition between qualita- 
tively different forms of relaxation behavior for K-SAT 
and K-COL. For small values of 7 or a, we find expo- 
nential relaxation to zero energy. For larger values of the 
parameter, the relaxation becomes algebraic. For large 
enough 7 or a, the energy freezes at a non-zero value. 

The change between exponential and power law relax- 
ation occurs for smaller values of 7 and a than the tran- 
sition in solvability. More extensive numerical investiga- 
tions will be necessary to determine the exact nature of 
this transition, in particular whether it is a sharp transi- 
tion, or whether there are intermediate regio ns wi th dif- 
ferent forms of relaxation behavior. In section VII below, 



FIG. 1. This graph is 3-colorable, but the Monte Carlo 
algorithm can get stuck in the local minimum shown (letters 
denote colors). 



we study the freezing transition, where the final energy of 
the MC algorithm changes from zero to a non-zero value. 
When this quantity is measured, a sharp transition is ob- 
served. This transition does not necessarily coincide with 
the change in functional form of the decay. 



In each time step in our simulations N spin flips are 
attempted. Each flip consists of selecting a random spin 
and randomly changing its color. The flip is allowed if 
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t (MC steps per spin) 

FIG. 2. The energy per spin as a function of time averaged 
over two graphs of size 10 7 for 3-COL. From bottom left to 
top right, 7 = 1,2, 3, 4, 5, 6, 7, 8. 



t (MC steps per spin) 



FIG. 3. The energy per spin as a function of time for 
3-SAT with 10 6 variables, averaged over 10 different formu- 
las. From bottom left to top right, a = 2, 3, 4, 6. 

Figure || shows the relaxation behavior for 3-COL in 
a log-log diagram, while figure |^ shows that the behav- 
ior for 3-SAT is similar. The data in these figures was 
determined for systems of size 10 7 (for 3-COL) and 10 6 
(for 3-SAT). The relaxation showed similar behavior for 
other values of K as well. 



The change between exponential and power law decay 
is illustrated in more detail for 3-COL in figure |], where 
we have plotted e(t) — e(500) for systems of N = 10 6 
spins. The figure shows data for 7 = 0.5 up to 3.0 in 
increments of 0.1. For small 7, the decay is exponential 
(see also figure || below for 7 = 1). When 7 is increased, 
a change from exponential to power law behavior is ob- 
served. A reasonable fit to power law behavior is found 
approximately for 7 > 2. 
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FIG. 5. Connectivity 7 = 1 (lower curve) in 3-COL gives 
a reasonable fit to an exponential exp (— 1/2) for 5 < t < 22. 
Data for 7 = 1.5 (upper curve) is also shown. 

Figure |^ shows the data for 7 = 1 , where a reasonable 
fit to exponential relaxation is obtained, and for 7 = 
1.5, where a cross-over behavior is seen (these and the 
following data were obtained for systems of size N = 10 7 ) 
Figure || shows the data for 7 = 2, 3, and 4, where a 
power law e ~ eo + t 11 is found. Figure Q shows that 
the power law also applies for 7 = 5 and 8; here the 
exponent is given by /i « 0.85. 





1 (MC steps per spin) 

FIG. 4. The energy per spin as a function of time, with 
the energy at 500 MC steps subtracted, for 3-COL with 
7 = 0.5, 0.6, . . . 3.0. The figure illustrates the transition from 
fast to power law relaxation. 



t (MC steps per spin) 



FIG. 6. For 3-COL, subtracting a constant from the en- 
ergy gives power law relaxation for 7 = 2 (bottom curve), 3 
(middle curve) and 4 (top curve). Note the finite-size effects 
in the bottom curve. 
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FIG. 7. For 3-COL, subtracting a constant from the en- 
ergy gives power law relaxation with approximately identical 
exponents n m 0.85 for 7 = 5, 6, 7 and 8. Data is plotted for 
7 = 5 (lower curve) and 7 = 8 (upper curve). 

The exponents for 3-COL are summarized in table Q; 
note that due to finite-size effects only data up to t = 200 
was used to determine the exponent for 7 = 2. 

The values for which the energy freezes for different 7 
in 3-COL are shown in figure ^|. For large 7, an approx- 
imately linear increase is observed. 

For 3-SAT, a — 2 gives exponentially fast decay, while 
there appears to be a crossover behavior for a = 3. For 
a = 4 and 6, power law relaxation is obtained, e ~ e + 
i~ M with \i k, 0.6 in both cases, see figure ^|. As in 3- 
COL, eo was found to increase approximately linearly 
with a. 

No significant change in the behavior was seen in finite 
temperature simulations. Raising the temperature makes 
it possible to escape from one local minimum, but the 
system is then trapped in another before the ground state 
is reached. Raising the temperature further repeats this 
scenario but also increases the fluctuations in the energy. 
For high enough temperatures, the fluctuations take over 
completely and the system is not trapped in any local 
minimum. 
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TABLE I. Approximate exponents for 3-COL. The relax- 
ation behavior can be described by e ~ eo + i _A \ with different 
constants eo for different 7's. The exponents were determined 
using two graphs of size 10 7 . 



FIG. 8. The value at which the energy freezes as a function 
of 7 for 3-COL with 10 7 variables. 
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FIG. 9. For 3-SAT, power law relaxation is obtained for 
a — 4 (bottom curve) and 6 (upper curve). The same data 
is shown as in figure ^, but with an a-dependent constant 
subtracted from the energy. 

We found the same form of power law relaxation with 
approximately identical exponents for temperatures up 
to T — 0.4, but the frozen- in value of the energy, eo, 
depended on the temperature, see figure [l^. For 7 = 4, 
the T = 0.2 runs were able to achieve an almost 30% 
better state than the T = runs. 



FIG. 10. The energy after 10 3 MC updates per spin as a 
function of temperature for 7 = 2 to 7. 



Using simulated annealing to find the ground states 
of some NP-complete spin glasses (3D ± J and infinite 
range models), Grest et al [Q found a logarithmic decay 
of the energy, e ~ eo + 1/logi. Similar decay was also 
recently found by Kiihn et al jl5| , who used an algorithm 
where attempts were made to flip several spins at once. 
The number of simultaneous flips, which plays the same 
role as the temperature in simulated annealing, was then 
slowly decreased. These methods will always find the 
ground state, whereas the T = MC algorithm studied 
here can get stuck in local minima. The faster relaxation 
of the T — MC method thus comes at the price of 
having no guarantee of finding the ground state. 

The exact values of the critical parameter for K-COL 
varies depending on the ensemble of graphs used || . We 
tried different ensembles and found no significant differ- 
ences in the relaxation behavior (e.g., 0™-model Feynman 
diagrams show approximately the same behavior as ran- 
dom graphs with connectivity 7 = n). 

Since the MC algorithm can get stuck in local minima, 
one should verify that the behavior does not depend on 
the initial values of the spins. It is also important to 
check to what extent the result depends on the choice of 
random graph. To test this, we have performed simula- 
tions where in addition to averaging over several different 
graphs we also restarted the MC algorithm with differ- 
ent inital spin configurations. In particular, 3-COL with 
7 = 2 and 4 was studied in detail. 




FIG. 11. For 3-COL the mean energies of several runs with 
(N a ,N r ) = (1,1000) are compared to the results for N = 10 7 
from above for 7 = 2 (lower curves, 56 different runs) and 
7 = 4 (upper curves, 54 different runs). 

In figure [ll] we compare the results for N = 10 7 shown 
above to runs for a smaller system (N = 10 4 ), where an 
average over 1000 initial states was performed for each 
of a larger set of graphs. The figure shows the average 
energy from these runs and e(t) from the runs with N = 
10 7 for 7 = 2 and 4. A reasonable agreement is found in 
both cases. 

The variation of the result depending on the choice of 
graph and initial state is further illustrated in figures |lj 
and [l3|, which show the average energy and the standard 
error 



(where (■) denotes the average over graphs and restarts, 
and N g and N r stand for the number of different graphs 
and the number of runs, respectively) for N — 10 4 and 
7 = 2, for each of 56 different runs with N r = 1 and N g 
1000. Figure [l4| shows the energy and standard error for 
N = 10 4 and 7 = 4 for 54 different runs with N r = 1 
and N g = 1000, while figure |l5| shows the standard errors 
from each of these runs. 
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FIG. 12. The energy as a function of time for 56 runs with 
(N g ,N T ) = (1,1000), for 3-COL with 7 = 2. 
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FIG. 13. The standard error as a function of time for 56 
runs with (N g ,N r ) = (1,1000), for 3-COL with 7 = 2. 

For 7 = 4, the variation among runs on a single graph 
(see figure Eq) is very small compared to the variation 
among randomly generated graphs shown in figure [l4|; 
for 7 = 2 these are of comparable magnitude. 
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t (MC steps per spin) 

FIG. 14. The energy as a function of time for 54 runs with 
{N g ,N r ) = (1,1000), for 3-COL with 7 = 4. 




t (MC steps per spin) 

FIG. 15. The standard error as a function of time for 54 
runs with {N g ,N r ) = (1,1000), for 3-COL with 7 = 4. 

We also found similar relaxation behavior when we 
started with all spins having the same value, and did 
not find any differences using different random number 
generators. The generators used include the Mitchell- 
Moore additive generator (e.g., p6[), a 48-bit multiplica- 
tive congruential generator with multiplicative and addi- 
tive constants 25214903917 and 11, respectively, and the 
standard C library randO function. 



VII. COMPARISON TO OTHER LOCAL SEARCH 
METHODS 

The Monte Carlo method is an example of a local 
search method. Local search methods start with a candi- 
date solution and try to improve it by changing it locally 
(e.g., by flipping a spin). Various other such methods 
have been used to study the phase transition in search 
cost in NPC-problems. Gent and Walsh Q studied 
small problems using the GSAT method, a hill-climbing 
procedure for solving K-SAT problems, and found ex- 
ponential relaxation over very small time ranges. The 
GSAT method is similar to MC simulations, but differs 
in the way a spin is chosen to flip. Instead of flipping 
a random spin, GSAT selects that spin whose flipping 
will decrease the energy the most. Clark ct al Egf used 



only solvable problems (determined by first using a com- 
plete backtracking method) and found an easy-hard-easy 
transition in search cost for two local search methods. 
The hardest problems occurred approximately at the true 
transition. 




(MC steps per spin) 



FIG. 16. The relaxation in 3-COL with only colorable 
graphs allowed. No significant change compared to figure || is 
observed. From bottom to top, 7 = 1,2, 3, 4, 5, 6, 7 and 8. 

In our case, the behavior is quite different. This is par- 
ticularly evident in figure 16, which shows the relaxation 
for 3-COL with 10 5 spins using only colorable graphs. 
It shows the same behavior as in figure |^ for both col- 
orable and uncolorable problems. An exceptional case 
is 2- COL, where we found different behavior when only 
colorable graphs were used, see figure [l7]. The reason for 
this is that for K = 2 and colorable graphs there is no dif- 
ference between the ferromagnetic and antiferromagnetic 
models. Results and arguments for the relaxation behav- 
ior of ferromagnetic Potts models on random graphs will 
be presented elsewhere pOA. 



t (MC steps per spin) 



FIG. 17. The relaxation in 2-COL with only colorable 
graphs allowed. The behavior is similar to that found for 
ferromagnetic random graphs. 7=1 starts with the lowest 
energy, followed by 7 = 2,3,4,5,6,7,8. The small arrow indi- 
cates the 7 = 8 data. 

Problems with large connectivities are thus always 
hard to solve using the Monte Carlo method. This is a 
short-coming of the MC algorithm — even if the problem 
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is solvable, the MC algorithm can get stuck in local min- 
ima that the other, smarter local search methods manage 
to avoid. 

In order to quantify the difference between the MC al- 
gorithm and other local search methods, and to compare 
the freezing transition with the solvability transition, we 
determined the fraction rj(t) of problems for which the 
MC method did not find an e = ground state in less 
than t MC steps per spin. We found that r\ displayed a 
behavior similar to the fraction of solvable problems - 
there appears to be a phase transition, but for a smaller 
value of 7. 




-2-1 1 2 3 4 5 



FIG. 18. Finite-size scaling analysis of the fraction of prob- 
lems for 3-COL where a zero energy ground state has not 
been found before 10 5 MC updates per spin. The system size 
N ranges from 20 to 1000, and the data is plotted against 
the rescaled parameter (7/7,3 — lJAT 1 /", with j c = 2.4 and 
v = 3.75. A good fit is obtained, except for N = 20 and 30. 

In figure U, we plot r){10 5 ) for 3-COL against the 
rescaled parameter 

- l)N^ v (7) 

with 7 C = 2.4 and v — 3.75 and N ranging from 20 to 
1000. It is clear that there is a freezing transition well be- 
low the occurrence of the solvability transition. However, 
it is necessary to be cautious when drawing conclusions 
from small systems. The largest systems we have simu- 
lated are of size 10 7 for 3-COL and 10 6 for 3-SAT; these 
systems were used to fit the relaxation behavior shown in 
table I above. To determine the freezing transition, we 
used considerably smaller systems. 

The approximate values where the transition in fig- 
ure |l8| was found are shown in table |l| for 2 < K < 5, 
together with the experimental values for the solvability 
transition for K-SAT from @, and values for the solv- 
ability transition for K-COL. The latter were obtained 
by ourselves using a non-optimized backtrack-search pro- 
gram with the Brelaz heuristic j49|. For K — 3, our value 
coincides with the literature [|| ; we are not aware of any 
values for K = 4 and 5 in the literature. For K > 4, 
the values for the solvability transition are probably not 



very accurate — the results vary strongly with the num- 
ber of variables and the number of graphs tested. It is 
clear, however, that the relaxation transition happens for 
smaller values than the solvability transition. 

Most of our data for the freezing transition was de- 
termined by averaging over between 100 (for N = 1000) 
and 1000 different graphs with a single MC run on each 
graph. We have also made some runs where we restarted 
the MC algorithm using different initial spin configura- 
tions for each graph. For small N, rj(t) varied about 13% 
depending on whether we used 1, 10 or 100 restarts, but 
for N > 80 there was essentially no difference between 
the different runs, see figure |l^. In each of these runs we 
averaged over 1000 different graphs. 



0.6 - I 
1(1° 4 ) 0.5 - 



FIG. 19. The fraction 77 of problems for 3-COL where a 
zero energy ground state has not been found after 10 4 MC 
updates per spin, averaged over 1000 different graphs of size 
100, with 1, 10 or 100 restarts per graph. 

The precise location of the freezing transition depends 
on the algorithm used to determine T)(t). For example, 
finite temperature MC methods give a higher value for 
the location of the transition. This illustrates the differ- 
ence between the freezing transition and the solvability 
transition. The former is a property of the method used 
to solve a problem, while the latter is a characteristic of 
the problem itself. 



K = 


2 


3 


4 


5 


K-COL 

Ic 


0.2 


2.4 
4.6 


6.7 
8.7 


11 
13.1 


K-SAT 

a e c xp 


1 

1.0 


4 
4.17 


8.4 
9.75 


16.7 
20.9 



TABLE II. The approximate values of 7 and a where the 
freezing transition was found, determined using a simple FSS 
scaling. Also shown are experimental values y^ mp and a™ 11 
for the solvability transition determined by Kirkpatrick and 
Selman and by us. 
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VIII. OTHER RESULTS 

We have also found interesting scaling relations for 
other quantities. For graph coloring, the energy as we 
have defined it above measures the fraction of frustrated 
edges. Another possibility is to measure the fraction of 
frustrated spins, e s . This quantity shows the same re- 
laxation behavior, if used in the Metropolis algorithm 
instead of e. 

It is also interesting to study the ratio 7f rus t = — , 
which is the connectivity of the sub-graph spanned by 
the frustrated edges. Figure [2C| suggests a transition from 
7frust = 1 for easy problems to 7f rus t > 1 for harder. If the 
connectivity is 1, all frustrated edges are isolated, while 
connectivity 2 would mean that all frustrated spins were 
connected in chains. 




FIG. 20. The connectivity of the graph spanned by the 
frustrated edges as a function of time. Note the change from 
7 frust = 1 to > 1 as the problem difficulty increases. The data 
are for 3-COL with N — 10° , averaged over 10 graphs. From 
bottom to top, 7 = 2, 3, 4, 5, 6, 7. 

The Monte Carlo dynamics itself has interesting prop- 
erties. Figure |l] shows the fraction r(t) of persistent 
spins, i.e., those that have not yet been flipped, as a 
function of time for 3-SAT. The data suggest a transi- 
tion from an exponentially fast to a logarithmically slow 
decay as a is increased (this will be explored further else- 
where). For the Ising and Potts models on a square lat- 
tice, this quantity has been found to scale with a power 
law dependence on time [pi 52[. 




(MC steps pet spin) 



FIG. 21. For 3-SAT, the fraction of persistent spins shows 
a similar transition as the energy. The data are averages of 
10 runs on systems of 10 6 spins and show a — 2 (bottom), 3, 
4, and 6 (top). 



IX. HEURISTIC ARGUMENTS FOR POWER 
LAW RELAXATION 

In K-COL and K-SAT there are global barriers to 
local improvement. In order to lower the frustration for 
a spin, we have to make changes to many other, unfrus- 
trated spins (compare figure pi). One way of explaining 
the relaxation is to look at these barriers and see how 
long it takes to overcome then. 

In the following, we call a proposed assignment of the 
N variables a state in the space of all solutions. Con- 
sider the probability p(i) that a state has i unsatisfied 
constraints, so that e = i/N. In analogy with expres- 
sion d) for K-COL, this can be approximated by 



M 

i. 



(8) 



where p v is the probability that a constraint is violated. 
The binomial factor ^ ^ represents the number of dif- 
ferent ways to choose the i unsatisfied edges. Recall that 
for K-COL p v = l/K and for K-SAT p v = 1/K 2 . In 
the approximation, we neglect all correlations between 
the constraints, such as the presence of triangles and 
other regular structures in the graph to be colored. In 
the following, only K-COL will be considered. 

The Monte Carlo algorithm works by generating a new 
state, and changing to this state if its energy is lower than 
or equal to that of the current state. In the approxima- 
tion we assume that the probability that the new state 
has energy e' is proportional to p S oi(e')- Given a state 
with energy e, all transitions to states with energy e' > e 
are forbidden. The probability of staying in a state with 

energy e is then proportional to YliLe P soi (*) > wnne t ne 
probability for a transition to an energy e" < e is pro- 
portional to p so \(e"). 

This process can be viewed as a Markov chain. The 
transition matrix between states of different energy has 
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elements pij given by the probability that state i is fol- 
lowed by state j: 

Pv = { E^HiPsol(n) = 1 - Eri=0ftol(n), i = j (9) 



Psol(j), 



I > J 



for i,j = 0...M. The transition probability pa can 
be written in terms of the hypergeometric function 
iT\ (a, b, c; z) |||| as 

Pit =Psoi(«)2^ r i(l,«- M, i+ 1; v ). 

Pu - 1 

This approximation ignores the detailed mechanisms of 
the MC-algorithm (which can only change a state lo- 
cally), instead we consider transitions between classes of 
states with the same energy. 



1= e 

2.8 ■ 
2.7 - 
2.6 ■ 
2.5 ■ 
2.J 
2 . 3 ■ 



FIG. 22. Log-log plot of the behavior of the Markov chain 
defined by equation ([]) . 

Using this transition matrix to evolve the states in time 
numerically, power law relaxation behavior is found, see 
figure |2^. The exponent differs from that in our MC 
simulations. One reason for this is that the MC sim- 
ulations only allow single-spin flips, while equation (j^) 
allows jumps between arbitrary spin configurations - 
the global and the local energy landscapes are different. 
We have performed simulations using a MC algorithm 
that changes the entire spin configuration instead of just 
a single spin. The results indicate a slower power law 
relaxation, with an exponent that is consistent with fig- 
ure ^2]. In contrast to the regular MC simulations, the 
exponent here showed a weak dependence on system size. 
The dependence was the same for both the simulations 
and time-evolution using equation (^|). 

In the MC algorithm, N attempts to flip a spin are 
made in each time step. Inspired by statistics of the 
energy landscape, we have made simulations where we 
require that there are N accepted spin flips in each time 
step (in T = simulations, this means that time is not 
increased for flip attempts to higher energy.) Time was 



increased by X/N for each change in the spin configura- 
tion. The power law relaxation behavior did not change 
when this modified algorithm was used. 

For these simulations, we can ignore all transitions to 
states with higher energy. Since a single spin flip is un- 
likely to decrease the energy by more than one (see fig- 
ure pp| ) , we can also ignore transitions to states with en- 
ergy j < i — 1, and consider a two-state system. 

To estimate the time needed to go from state i to j = 
i — 1, we assume that p(i)/p(i — 1) can be approximated 
by 



Psol(i)/Psol(« - 1) = 



p v M-i + 1 1 
1 — Pu i e 



(10) 



Since the average time needed to pass from state i to 
state j = i — 1 is given by the inverse of the transition 
probability per unit time, we then have 



At = 



p(i)+p(i - 1) 
p(i - 1) 



= 1 



Psoiji) 
Psol(* - 1) 



1+- (11) 



for the time needed to go from state i to i — 1. 

Another viewpoint is to consider the free energy barri- 
ers, A/, in the problem. Since transitions to states with 
higher energy are not allowed by the T — Monte Carlo 
dynamics, we can neglect the energy part of A/. The 
entropic part will however be positive, 



A/ - -k b T log N(j) + k b TlogN(i) = k b T log 



Psoiji) 

Psol(j') 

(12) 



for a transition from i to j. The time needed to overcome 
this free energy barrier is 



exp (— ) - Pso1 ^ ^ I 



(13) 



if i = j + 1 and e = i/N. 

By inverting equations (^lj) and ([31]), power law re- 
laxation with an asymptotic t~ x dependence can be ob- 
tained. This is justified for large 7, where all free energy 
barriers have roughly the same size, but not for 7 be- 
low the freezing transition. Similar difficulties also arise, 
e.g., for coarsening, where arguments for relaxation usu- 
ally are given for nucleation, where a single domain grows 
or shrinks. It is then justifiable to invert the equation re- 
lating the energy of the domain and the time needed to 
shrink it. In the coalescence regime, however, several 
such domains form and grow separately before they coa- 
lesce. The inversion is then not allowed, even though the 
same relaxation behavior is found in simulations. 

The power law with t~ l dependence approximately 
matches the behavior of K-COL for critically con- 
strained 7's (see table ||). For overconstrained problems, 
the relaxation is slower than f -1 . This could depend 
on the local energy landscape being different from the 
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global when the problems are overconstrained. The ar- 
guments above should be viewed only as a simplistic first 
approximation, and cannot be expected to reproduce all 
important features of the numerical results. 

The barriers in equation (|l2|), which are proportional 
to the logarithm of the energy, do not appear in any of the 
coarsening classes introduced by Lai et al j54|. The dif- 
ference seems to be that here there are entropy barriers. 
That the entropy is high for low energies is explained by 
the fact that the interactions are antiferromagnetic and 
the variables are Potts spins. If each spin can have K 
different values and the interactions are ferromagnetic, 
each satisfied bond can be satisfied in K different ways, 
while there are K 2 — K choices for the unsatisfied ones. 
If the interactions are antiferromagnetic, each satisfied 
bond has K 2 — K possibilities and the unsatisfied K. 
Since there are more satisfied bonds than unsatisfied, an 
antiferromagnetic problem is less restricted, at least for 
K > 2. As always, this reasoning ignores all correlations 
between constraints in the problem, such as triangles in 
graphs. 

X. CONCLUSIONS AND DISCUSSION 

We have shown that there are additional transitions in 
the behavior of two NPC problems — graph colorabil- 
ity and satisfiability. There is a freezing transition, and 
the relaxation changes from exponential to power law as 
the difficulty of the problem increases. These transitions 
occur for smaller values of the parameters 7 and a than 
the transitions in solvability and search cost. 

We also studied the effects of using only colorable 
graphs, and measured other quantities, such as the frac- 
tion of persistent spins for K-SAT, which also showed a 
transition. 

A simple heuristic argument for the form of the relax- 
ation behavior depending on entropy barriers was also 
given. The free energy barriers in these problems have a 
less pronounced energy dependence than in other mod- 
els. One can compare them to systems where logarithmic 
relaxation has been found for T ^ and T = MC sim- 
ulations, such as in fl55| , where a tiling model containing 
no randomness at all was studied. The main difference 
between these systems and ours seems to be that they 
have local interactions and frustration, whereas we have 
infinite-range interactions and frustration. 

The problems with logarithmic relaxation are in some 
ways simpler than the NPC problems studied here. 
Their ground states can always be found in polynomial 
time, and they do not contain any difficult optimization 
problem. On the other hand, K-SAT and K-COL have 
ground states that may require exponential time to find. 
Our results show that it can be easier to find approxi- 
mate solutions to these hard problems than for the easy 
problems giving logarithmic relaxation. Of course, to ac- 
tually find the ground state with the MC algorithm is 



very difficult for K-COL and K-SAT. What the NPC 
problems gain in finding an approximate solution is lost 
when it comes to finding the ideal, ground state solution. 

Other problems where entropy barriers and logarith- 
mic relaxation are found include the backgammon model 
introduced by Ritort |]56f and the random-field Ising 
model p7| , p8[ . The backgammon model contains no en- 
ergy barriers at all, only entropy barriers, and can be 
solved J59| by considering random-walk models with en- 
tropy barriers. In our models, the entropy barriers scale 
as the logarithm of the energy. This means that the time 
needed to overcome the free energy barrier loses its ex- 
ponential dependence and becomes algebraic, leading to 
power law relaxation. We have also found this behav- 
ior in more general constraint satisfaction problems than 
K-COL and K-SAT |5(|. Relaxation behavior similar 
to that in our models has also been obtained by Campel- 
lone et al for a short-range p-spin glass model |30| . These 
models should be studied further. 

We see several possibilities for future work in this field. 
More detailed studies of the behavior of the fraction of 
persistent spins as well as of damage spreading and hys- 
teresis in strongly disordered systems are needed. Since 
hard optimization problems are also spin glasses, a search 
for spin glass characteristics such as aging and chaotic 
responses to small perturbations in the hamiltonian may 
also be interesting. We are preparing a study of the en- 
ergy landscapes of these and related models. Investigat- 
ing the connection between the topology of the random 
graph and that of the energy landscape may also provide 
more insight into why NPC problems are difficult. 
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